function eq = solveEquilibrium(z0,lsap,pArrLength)

% Inputs:
% initial productivity z0
% size of the LSAP policy:  lsap = etag*(wg-1)
% lsap = etag*(wg-1); % Size of the lsap program

% The last argument is optional: Controls the length of pArr
% Default is 100. Set it to a higher level for improved resolution
if nargin<3
    pArrLength=100;
end

% Returns: An equilibrium structure with (among other things):
% The equilibrium price p and the interest rate rf (if multiple equilibria, return the best equilibrium)
% The required Sharpe ratio curve (for plotting purposes)
% The actual Sharpe ratio curve with rf=0 (for plotting purposes)


assert(lsap>=0);  % We are only interested in cases with positive lsap
assert(lsap<1); % The case above 1 doesn't make much sense
assert(z0<=1);  % We are only interested in cases with productivity shocks are below baseline
eq.lsap = lsap;
eq.z0 = z0;

% Underlying parameters(loaded through globals)
global g rho sigma PStar taub tauh k l tau tauBarBench

ph = exp(rho + g - sigma^2*(1-lsap)/tauh);
tau = taub*k+tauh*(1-k);
% Store all parameters for alter access
eq.params.g = g; eq.params.rho = rho; eq.params.sigma = sigma;
eq.params.tauh =tauh; eq.params.taub = taub;
eq.params.tau = tau;

eq.ph=ph;
% Baseline price without a demand recession
% Will adjust this for cases with demand recession
p = 1;

if ph >= 1
    % Simple case where households are sufficiently risk tolerant (with lsap)
    % No recession no matter what
    eq.unique=1;
    p=1;
    % No need to change the baseline levels
else
% ph < 1
% There can be demand recessions
zBar = fsolve(@(z) tau - k*l/z*(taub-tauh) -  sigma^2*(1-lsap)/(rho+g),1);
eq.zBar = zBar;

zTilde = l/ph;
eq.zTilde = zTilde;

% Set the appropriate arrays (for plotting purposes)
pArr = linspace(ph,1,pArrLength); % Normalized price array
eq.pArr = pArr;

sharpeActualArr = (g + rho - log(pArr))/sigma;
eq.sharpeActualArr = sharpeActualArr; 
% Calculate the Sharpe ratio also for the case when shock is temporary

% Calculate banks' implied balance sheets, and the implied risk tolerance, for each p
alphaArr = max(0,k*(1-l./(z0*pArr)));
eq.alphaArr = alphaArr;
tauBarArr = alphaArr*taub+(1-alphaArr)*tauh;
eq.tauBarArr = tauBarArr;

sharpeRequiredArr = sigma*(1-lsap)./tauBarArr;
eq.sharpeRequiredArr = sharpeRequiredArr;

% Check condition for uniquenes
if zTilde < zBar
   eq.unique = 1;
   if z0 <= zTilde
       % Corner solution where bank is bankrupt
       p = ph;
       ind = 1;
   elseif z0 >= zBar
       % Corner solution without demand recession
       p=1;
       ind = length(pArr);
       % 
   else
   % Find the unique intersection
   [sol,numCross] = findCross(sharpeRequiredArr,sharpeActualArr,pArr);
   assert(numCross==1); % Intersection is supposed to be unique
   p = sol{1}.x;
   ind = sol{1}.xind;
   % eq.Sharpe = sol.y;
   end
else
    % Case with possible multiplicity
    % Find all crossings (including corner solutions)
   [sol,numCross] = findCross(sharpeRequiredArr,sharpeActualArr,pArr);
   if numCross == 1
       eq.unique=1;
       % still unique
       p = sol{1}.x;
        ind = sol{1}.xind;
        % eq.Sharpe = sol.y;
   elseif numCross > 1
       % Multiple equilibria
       eq.unique=0;
       % Pick the best equilibrium (highest p)
       % Later, we can change to keep track of all equilibria. Not
       % necessary for now
       p = sol{numCross}.x;
       ind = sol{numCross}.xind;
   end
end
   
end

% Set the key equilibrium output
% The key thing that we should get right from the earlier analysis is p
% Given p (if correctly calculated), we can recover all else
alpha = max(0,k*(1-l/(p*z0)));
tauBar = alpha*taub+(1-alpha)*tauh;
rf = max(0,rho + g - (1-lsap)*sigma^2/tauBar); % Accommodates corner solutions
eq.p = p;
eq.ind = ind;
eq.alpha= alpha;
eq.tauBar = tauBar;
eq.rf = rf;


function [sol,numCross] = findCross(y1Arr,y2Arr,xArr)
    % To accommodate corner solutions, choose y1Arr and y2Arr so that
    % if y1Arr<y2Arr throughout then xArr(end) is the solution (and vice versa)
    
    % Find all crossings of two relations
    % y1(x1) and y2(x2)
    % defined on the same array, xArr
    % Requirement: The relations are continuous so every crossing will be interpreted as an equilibrium
    % Requirement: y1Arr(xArr)
    % Output:
    % sol{i}.y and sol{i}.x have the solutions corresponding to crossing i
    % sol{i}.xind has the x-index corresponding to the solution with crossing i
    % numCross: Number of equilibria      
    %
    % Corner cases: 
    %%%% If y1(1)>y2(1): Corner equilbirium at index 1
    %%%% If y1(end)<y2(end): Corner equilbrium at index end
    numCross =0; 
    ind1Arr = []; 
    lastDif = 0; % Keeps track of the last difference, y1-y2, at the last step
    % newDif will keep track of the difference at this iteration
    
    % Deal with for the low end explicitly to account for a corner solution
    newDif = y1Arr(1) - y2Arr(1);
    if newDif>=0
        display('Corner solution found!');
        numCross = numCross+1;
        sol{numCross}.y = y1Arr(1);
        sol{numCross}.x = xArr(1);
        sol{numCross}.xind = 1;
    end
    lastDif=newDif;
    % Loop for an interior solution
for ind=2:length(xArr)
    x = xArr(ind);
    y1 = y1Arr(ind);
    y2 = interp1(xArr,y2Arr,x,'linear');
    newDif = y1-y2;
    if isnan(y2)
        % No intersection
        continue;
    elseif newDif==0
        display('Crossing found!');
        numCross = numCross+1;
        sol{numCross}.y = y1;
        sol{numCross}.x = x;
        sol{numCross}.xind = ind;
    elseif newDif<0 & lastDif>0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        % Find the exact solution
        numCross = numCross+1;
        % Do a linear interpolation based on the magnitude of the differences
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    elseif newDif>0 & lastDif<0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        numCross = numCross+1;
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    end
    lastDif =newDif;
end  

    % Deal with for the high end explicitly to account for a corner solution
    N = length(xArr);
    newDif = y1Arr(N) - y2Arr(N);
    if newDif<=0
        display('Corner solution found!');
        numCross = numCross+1;
        sol{numCross}.y = y1Arr(N);
        sol{numCross}.x = xArr(N);
        sol{numCross}.xind = N;
    end

end

end